#coding=utf-8
import numpy as np
from numba import jit
from numba import cuda,float32,f4,void,int32
import time

def get_patch(input_array, i, j, kernel_width,kernel_height, stride):
    '''
    获取数组的切片
    '''
    #print input_array.shape
    if input_array.ndim==3:
        return input_array[:,i*stride:i*stride+kernel_width,j*stride:j*stride+kernel_height]
    if input_array.ndim==4:
        return input_array[:,:,i*stride:i*stride+kernel_width,j*stride:j*stride+kernel_height]
    else:
        return input_array[i*stride:i*stride+kernel_width,j*stride:j*stride+kernel_height]

def conv(input_array, kernel_array, output_array, stride, bias,axi):
    '''
    计算卷积，自动适配输入为2D和3D的情况
    '''
    output_width = output_array.shape[2]
    output_height = output_array.shape[1]
    kernel_width = kernel_array.shape[-1]
    kernel_height = kernel_array.shape[-2]
    for i in xrange(output_height):
        for j in xrange(output_width):
            output_array[:,i,j] = (get_patch(input_array, i, j, kernel_width,kernel_height, stride) * kernel_array).sum(axis=axi) + bias 
            
def padding(input_array, zp):
    '''
    为数组增加Zero padding，自动适配输入为2D和3D的情况
    '''
    if zp == 0:
        return input_array
    else:
        if input_array.ndim == 3:  # 三维数组
            input_width = input_array.shape[2]
            input_height = input_array.shape[1]
            input_depth = input_array.shape[0]
            padded_array = np.zeros((
                input_depth, 
                input_height + 2 * zp,
                input_width + 2 * zp)).astype('float32')
            padded_array[:,
                zp : zp + input_height,
                zp : zp + input_width] = input_array
            return padded_array
        if input_array.ndim == 4:  # 4维数组
            input_width = input_array.shape[3]
            input_height = input_array.shape[2]
            input_depth = input_array.shape[1]
            banch_size = input_array.shape[0]
            padded_array = np.zeros((
                banch_size,
                input_depth, 
                input_height + 2 * zp,
                input_width + 2 * zp)).astype('float32')
            padded_array[:,:,
                zp : zp + input_height,
                zp : zp + input_width] = input_array
            return padded_array        
        elif input_array.ndim == 2: #二维数组
            input_width = input_array.shape[1]
            input_height = input_array.shape[0]
            padded_array = np.zeros((
                input_height + 2 * zp,
                input_width + 2 * zp)).astype('float32')
            padded_array[zp : zp + input_height,
                zp : zp + input_width] = input_array
            return padded_array

# ################### GPU加速函数 #####################################

# 卷积核函数  （正向传播 ）
@cuda.jit(void(f4[:,:,:,:],f4[:,:,:,:],f4[:,:,:,:],f4[:],int32))
def fconvolution(arr,w,out,bias,stride):
    '''
    arr :输入的图片
    w: 卷积核
    out: 卷积后的图片
    '''
    filter_n,filter_c,filter_h,filter_w=w.shape  # 卷积核的数量，通道数 高 宽
    picture_n,picture_c,picture_h,picture_w=out.shape # 图片的数量，通道数 高 宽    
    x,y,z=cuda.grid(3) # x : 卷积后图片的第x 行 ；y: 卷积后图片的第y列; z: 第z张图片
    # 检查越界
    if x>=picture_h or y>=picture_w or z>=picture_n:
        return    
    # 载入图片数据
    #b=cuda.local.array((filter_c,filter_h,filter_w),dtype=float32) # 报错 不能动态声明数组
    b=arr[z,:,x*stride:x*stride+filter_h,y*stride:y*stride+filter_w]
    ws=w 
    bia=bias
    tmp=0.0
    for k in range(filter_n):
        tmp=0.0
        for i in range(filter_h):
            for j in range(filter_w):
                for c in range(filter_c):
                    tmp+=b[c,i,j]*ws[k,c,i,j]
                    #tmp+=arr[z,c,x*stride+i,y*stride+j]*w[k,c,i,j]
        out[z,k,x,y]=tmp+bia[k]

# 卷积核函数  （反向传播 ）
@cuda.jit(void(f4[:,:,:,:],f4[:,:,:,:],f4[:,:,:,:]))
def bconvolution(arr,w,out):
    '''
    arr :输入的图片
    w: 卷积核
    out: 卷积后的图片
    '''
    stride=1
    filter_n,filter_c,filter_h,filter_w=w.shape  # 卷积核的数量，通道数 高 宽
    picture_n,picture_c,picture_h,picture_w=out.shape # 图片的数量，通道数 高 宽    
    x,y,z=cuda.grid(3) # x : 卷积后图片的第x 行 ；y: 卷积后图片的第y列; z: 第z张图片
    # 检查越界
    if x>=picture_h or y>=picture_w or z>=picture_n:
        return        
    # 载入图片数据
    #b=cuda.local.array((filter_c,filter_h,filter_w),dtype=float32) # 报错 不能动态声明数组
    b=arr[z,:,x*stride:x*stride+filter_h,y*stride:y*stride+filter_w]
    ws=w 
    tmp=0.0
    for k in range(filter_c):  # 输出的图片有filter_c个通道
        tmp=0.0
        for i in range(filter_h):
            for j in range(filter_w):
                for c in range(filter_n):
                    tmp+=b[c,i,j]*ws[c,k,i,j]
                    #tmp+=arr[z,c,x*stride+i,y*stride+j]*w[k,c,i,j]
        out[z,k,x,y]=tmp

# max polling  （正向传播 ）
@cuda.jit(void(f4[:,:,:,:],f4[:,:,:,:],int32[:,:,:,:,:],int32,int32,int32,int32))
def maxpolling(arr,out,max_idx,filter_h,filter_w,stride_h,stride_w):
    '''
    arr :输入的图片
    out: 卷积后的图片   （输出）
    max_idx: 最大元素的索引数组（输出）
    filter_h: 窗口的高
    filter_w：窗口的宽
    stride_h,stride_w :两个方向的步长（一般等于核对尺寸）
    '''
    #filter_n,filter_c,filter_h,filter_w=w.shape  # 卷积核的数量，通道数 高 宽
    picture_n,picture_c,picture_h,picture_w=out.shape # 图片的数量，通道数 高 宽    
    x,y,z=cuda.grid(3) # x : 卷积后图片的第x 行 ；y: 卷积后图片的第y列; z: 第z张图片
    # 检查越界
    if x>=picture_h or y>=picture_w or z>=picture_n:
        return    
    # 载入图片数据
    #b=cuda.local.array((filter_c,filter_h,filter_w),dtype=float32) # 报错 不能动态声明数组
    b=arr[z,:,x*stride_h:x*stride_h+filter_h,y*stride_w:y*stride_w+filter_w]
    # 找出每张图片每个通道的最大值
    tmp=0.0
    for k in range(picture_c):
        tmp=b[k,0,0]
        ix=0
        iy=0
        for i in range(filter_h):
            for j in range(filter_w):
                    if tmp<b[k,i,j]:
                        tmp=b[k,i,j]
                        ix=i
                        iy=j
        out[z,k,x,y]=tmp #保存最大值
        max_idx[z,k,x,y,0]=x*stride_h+ix #保存最大值索引
        max_idx[z,k,x,y,1]=y*stride_w+iy #保存最大值索引

# max polling  （反向传播 ）
@cuda.jit(void(f4[:,:,:,:],f4[:,:,:,:],int32[:,:,:,:,:],int32,int32,int32,int32))
def bmaxpolling(arr,out,max_idx,filter_h,filter_w,stride_h,stride_w):
    '''
    arr :输入的图片
    out: 卷积后的图片   （输出）
    max_idx: 最大元素的索引数组（输入）
    filter_h: 窗口的高
    filter_w：窗口的宽
    stride_h,stride_w :两个方向的步长（一般等于核对尺寸）
    '''
    #filter_n,filter_c,filter_h,filter_w=w.shape  # 卷积核的数量，通道数 高 宽
    picture_n,picture_c,picture_h,picture_w=arr.shape # 图片的数量，通道数 高 宽    
    x,y,z=cuda.grid(3) # x : 卷积后图片的第x 行 ；y: 卷积后图片的第y列; z: 第z张图片
    # 检查越界
    if x>=picture_h or y>=picture_w or z>=picture_n:
        return    
    for k in range(picture_c):
        ix=max_idx[z,k,x,y,0]
        iy=max_idx[z,k,x,y,1]
        out[z,k,ix,iy]=arr[z,k,x,y]
        
# #############################################################

def fconv(input_array, kernel_array,bias,stride):
    '''
    input_array: 被卷积的图片 
    kernel_array： 卷积核
    output_array：输出
    stride：卷积步长
    bias：卷积偏置
    '''
    filter_n,filter_c,filter_h,filter_w=kernel_array.shape  # 卷积核的数量，通道数 高 宽
    picture_n,picture_c,picture_h,picture_w=input_array.shape # 图片的数量，通道数 高 宽
    
    # 计算卷积后图片的尺寸
    #stride=1 #卷积步长
    height=(picture_h-filter_h)/stride +1  # 应该确保能整除
    width=(picture_w-filter_w)/stride +1
    
    # 在显存建立一个空数组存储卷积后的图片
    out=cuda.device_array((picture_n,filter_n,height,width),dtype=np.float32)
    
    # 计算 block 和 gird 的尺寸 每个block有128个线程(每个block处理8张图片的一部分)
    block_x=4
    block_y=4
    block_z=8
    block=(block_x,block_y,block_z)
    
    grid_x=(height+block_x-1)/block_x
    grid_y=(width+block_y-1)/block_y
    grid_z=(picture_n+block_z-1)/block_z
    grid=(grid_x,grid_y,grid_z)
    
    da = cuda.to_device(input_array)  # 显式的将数组拷贝到显存能大大减小时间（不用拷贝回来）
    dw=cuda.to_device(kernel_array)    
    db=cuda.to_device(bias)
    fconvolution[grid,block](da,dw,out,db,stride)  # a,w隐式的拷贝到显卡的global memory
    
    # 从显存拷贝数组到主机内存
    output=out.copy_to_host()
    del out,da,dw,db # 释放显存
    return output

def bconv(input_array, kernel_array):
    '''
    input_array: 被卷积的图片(误差) 
    kernel_array： 卷积核
    output_array：输出
    stride：卷积步长
    bias：卷积偏置
    '''
    filter_n,filter_c,filter_h,filter_w=kernel_array.shape  # 卷积核的数量，通道数 高 宽
    picture_n,picture_c,picture_h,picture_w=input_array.shape # 图片的数量，通道数 高 宽
    stride=1
    # 计算卷积后图片的尺寸
    #stride=1 #卷积步长
    height=(picture_h-filter_h)/stride +1  # 应该确保能整除
    width=(picture_w-filter_w)/stride +1
    
    # 在显存建立一个空数组存储卷积后的图片
    out=cuda.device_array((picture_n,filter_c,height,width),dtype=np.float32)
    
    # 计算 block 和 gird 的尺寸 每个block有128个线程(每个block处理8张图片的一部分)
    block_x=4
    block_y=4
    block_z=8
    block=(block_x,block_y,block_z)
    
    grid_x=(height+block_x-1)/block_x
    grid_y=(width+block_y-1)/block_y
    grid_z=(picture_n+block_z-1)/block_z
    grid=(grid_x,grid_y,grid_z)
        
    #启动内核
    da = cuda.to_device(input_array)  # 显式的将数组拷贝到显存能大大减小时间（不用拷贝回来）
    dw=cuda.to_device(kernel_array)    
    
    bconvolution[grid,block](da,dw,out)  # a,w隐式的拷贝到显卡的global memory
    
    # 从显存拷贝数组到主机内存
    output=out.copy_to_host()
    del out,da,dw # 释放显存
    return output

def fmaxpoll(input_array,filter_h,filter_w,stride_h,stride_w):
    '''
    input_array: 被卷积的图片
    filter_h,filter_w,stride_h,stride_w: 卷积核相关参数
    '''
    picture_n,picture_c,picture_h,picture_w=input_array.shape # 图片的数量，通道数 高 宽
    # 计算卷积后图片的尺寸
    #stride=1 #卷积步长
    height=(picture_h-filter_h)/stride_h +1  # 应该确保能整除
    width=(picture_w-filter_w)/stride_w +1
    # 在显存建立一个空数组存储卷积后的图片
    out=cuda.device_array((picture_n,picture_c,height,width),dtype=np.float32)
    max_idx=cuda.device_array((picture_n,picture_c,height,width,2),dtype=np.int32)
    # 计算 block 和 gird 的尺寸 每个block有128个线程(每个block处理8张图片的一部分)
    block_x=4
    block_y=4
    block_z=8
    block=(block_x,block_y,block_z)
    grid_x=(height+block_x-1)/block_x
    grid_y=(width+block_y-1)/block_y
    grid_z=(picture_n+block_z-1)/block_z
    grid=(grid_x,grid_y,grid_z)
    #启动内核
    da = cuda.to_device(input_array)  # 显式的将数组拷贝到显存能大大减小时间（不用拷贝回来) 
    maxpolling[grid,block](da,out,max_idx,filter_h,filter_w,stride_h,stride_w)
    # 从显存拷贝数组到主机内存
    output=out.copy_to_host()
    outidx=max_idx.copy_to_host()
    del out,da,max_idx # 释放显存
    return output,outidx

def bmaxpoll(input_array,max_idx,filter_h,filter_w,stride_h,stride_w):
    '''
    input_array: 被卷积的图片
    filter_h,filter_w,stride_h,stride_w: 卷积核相关参数
    max_idx: 最大值索引数组
    '''
    #filter_n,filter_c,filter_h,filter_w=kernel_array.shape  # 卷积核的数量，通道数 高 宽
    picture_n,picture_c,picture_h,picture_w=input_array.shape # 图片的数量，通道数 高 宽
    
    # 计算卷积后图片的尺寸
    #stride=1 #卷积步长
    #height=(picture_h-filter_h)/stride_h +1  # 应该确保能整除
    #width=(picture_w-filter_w)/stride_w +1
    height=picture_h*filter_h  # 应该确保能整除
    width=picture_w*filter_w
    
    # 在显存建立一个空数组存储卷积后的图片
     # 注意！这样申请的数组不会初始化为0，如果在内核中没有对out的每个元素赋值，则结果是不可预测的！ 而这里的out恰好是这种情况，所以不能这样申请数组
    #out=cuda.device_array((picture_n,picture_c,height,width),dtype=np.float32)  # 程序bug的原因！
    out=np.zeros((picture_n,picture_c,height,width),dtype=np.float32)          # 这才是正确的（用host数组初始化为0）
   
    # 计算 block 和 gird 的尺寸 每个block有128个线程(每个block处理8张图片的一部分)
    block_x=4
    block_y=4
    block_z=8
    block=(block_x,block_y,block_z)
    
    grid_x=(picture_h+block_x-1)/block_x
    grid_y=(picture_w+block_y-1)/block_y
    grid_z=(picture_n+block_z-1)/block_z
    grid=(grid_x,grid_y,grid_z)
     
    #for i in range(1):
    #启动内核
    #da = cuda.to_device(np.ascontiguousarray(input_array))  # 显式的将数组拷贝到显存能大大减小时间（不用拷贝回来) 
    da = cuda.to_device(input_array)
    dmax=cuda.to_device(max_idx)
    bmaxpolling[grid,block](da,out,dmax,filter_h,filter_w,stride_h,stride_w)
    # 从显存拷贝数组到主机内存
    del da,dmax # 释放显存
    return out
# #######################################################################

# 卷积层类
class ConvLayer(object):
    def __init__(self,
                 filter_width, 
                 filter_height, filter_number, 
                 zero_padding, stride, activator,
                 learning_rate,input_shape):
        '''
        初始化卷积层
        '''
        self.filter_width = filter_width  # 卷积核的宽
        self.filter_height = filter_height # 卷积核的高
        self.filter_number = filter_number  # 卷积核的个数
        self.zero_padding = zero_padding   # 图片边缘填充0的个数
        self.stride = stride  # 卷积核移动的步长
        self.input_shape=input_shape
        self.channel_number,self.input_width,self.input_height=input_shape
        # 计算输出图片的宽
        self.output_width = (self.input_width- filter_width+2*zero_padding)/stride+1
        # 计算输出图片的高
        self.output_height=(self.input_height- filter_height+2*zero_padding)/stride+1
        # 储存输出图片
        #self.output= np.zeros((self.filter_number,self.output_height, self.output_width))
        # 初始化卷积核
        np.random.seed(0)
        self.weights = np.random.uniform(0, 1e-1,(self.filter_number, self.channel_number,self.filter_height, self.filter_width)).astype('float32') #一个四维数组，矢量化编程
        self.weights_bias=np.zeros(self.filter_number).astype('float32')  # filter_number个偏置
        self.weights_grad=np.zeros(self.weights.shape).astype('float32') # 权重的导数
        self.activator = activator
        self.learning_rate = learning_rate
        self.initialized=0
    def forward(self, input_array):
        '''
        计算卷积层的输出
        输出结果保存在self.output_array
        '''
        banch_size,channel_num,height,width=input_array.shape
        self.banch_size=banch_size
        # 如果没有初始化就进行初始化
        if not self.initialized:
            self.initialized=1
        self.input_array = input_array
        # 图片边缘补0
        self.padded_input_array = padding(input_array,self.zero_padding)   # 可以实现四维矢量操作
        self.output=fconv(self.padded_input_array.astype('float32'), self.weights.astype('float32'), self.weights_bias.astype('float32'),self.stride)
        self.output=self.activator.forward(self.output.astype('float32'))
        
    def create_delta_array(self):
        '''
        保存上一层的敏感度数组
        '''
        return np.zeros((self.banch_size,self.channel_number,
            self.input_height, self.input_width)).astype('float32')    
    
    def backward(self, sensitivity_array):
        '''
        计算传递到上一层的sensitivity map ( 用于计算上一层的权重梯度)
        sensitivity_array: 本层的sensitivity map 
        activator: 上一层的激活函数  
        '''
        # 乘以激活函数的倒数，获得本层的敏感度（本层神经元输入对误差的导数）
        sensitivity_array=self.activator.backward(self.output)*sensitivity_array  # 有问题?
        # 处理卷积步长，对原始sensitivity map进行扩展
        expanded_array = self.expand_sensitivity_map(sensitivity_array)
        expanded_width = expanded_array.shape[-1]
        #zp = (self.input_width + self.filter_width - 1 - expanded_width) / 2 # 无法理解 ，似乎是错的
        zp = self.filter_width - 1  # 这样似乎才是对的
        padded_array = padding(expanded_array, zp)
        #print 'padded_array',padded_array.shape
        # 边缘添加的0不计算敏感度，因此可以截取敏感度矩阵的一部分
        pw,ph=padded_array.shape[-1],padded_array.shape[-2]
        padded_array=padded_array[:,:,self.zero_padding:ph-self.zero_padding,self.zero_padding:pw-self.zero_padding]
        # 对于具有多个filter的卷积层来说，最终传递到上一层的
        # sensitivity map相当于所有的filter的sensitivity map之和
        # 将filter权重翻转180度
        flipped_weights =np.rot90(self.weights, 2, (2,3))
        self.delta=bconv(padded_array.astype('float32'), flipped_weights.astype('float32'))
        # 最终要传递到上一层的敏感度
        # 接下来 ，计算本层权值的梯度 expanded_array
        for i in range(self.filter_number):
            conv(self.padded_input_array, expanded_array[:,i:i+1,:,:],self.weights_grad[i], 1, 0,axi=(0,2,3))
        self.bias_grad = expanded_array.sum(axis=(0,2,3)) 
    
    def expand_sensitivity_map(self, sensitivity_array):
        depth = sensitivity_array.shape[1]
        # 确定扩展后sensitivity map的大小
        # 计算stride为1时sensitivity map的大小
        expanded_width = (self.input_width - 
            self.filter_width + 2 * self.zero_padding + 1)
        expanded_height = (self.input_height - 
            self.filter_height + 2 * self.zero_padding + 1)
        # 构建新的sensitivity_map
        expand_array = np.zeros((self.banch_size,depth, expanded_height, expanded_width)).astype('float32')
        # 从原始sensitivity map拷贝误差值
        for i in range(self.output_height):
            for j in range(self.output_width):
                i_pos = i * self.stride
                j_pos = j * self.stride
                expand_array[:,:,i_pos,j_pos] = \
                    sensitivity_array[:,:,i,j]
        return expand_array        

    def update(self,rate=0.1):
        '''
        按照梯度下降，更新权重
        '''
        self.weights_bias-=rate*self.bias_grad
        self.weights-=rate*self.weights_grad
            
class MaxPoolingLayer(object):
    def __init__(self,filter_width, 
                 filter_height, stride):
        self.filter_width = filter_width  #这里的filter不是卷积中的filter，而只是一个取样点区域。。。其实也差不多
        self.filter_height = filter_height
        self.stride = stride
        # 只是为了保持层的结构一致，没有实际用处
        self.weights=np.zeros(1)
        self.weights_bias=np.zeros(1)
        self.initializied=False
        self.input_shape=[]
        
    def forward(self, input_array):  # max polling    
        self.output,self.max_idx=fmaxpoll(input_array.astype('float32'),self.filter_height,self.filter_width,self.stride,self.stride)
        
    def backward(self,sensitivity_array):
        self.delta=bmaxpoll(sensitivity_array.astype('float32'),self.max_idx.astype('int32'),self.filter_height,self.filter_width,self.stride,self.stride)
                    
    def update(self, learning_rate):
        pass 

# 压平层
class FlattenLayer(object):
    def __init__(self):
        self.input_shape=None
        self.output=None
        self.delta=None
        # 只是为了保持层的结构一致，没有实际用处
        self.weights=np.zeros(1)
        self.weights_bias=np.zeros(1)
        self.initializied=False
    def forward(self,input_arr):
        self.input_shape=input_arr.shape
        self.output=(np.array(map(lambda a:a.flatten(),input_arr))).T
    def backward(self,delta_array):
        self.delta=np.ascontiguousarray((delta_array.T).reshape(self.input_shape),dtype=np.float32) #确保数组是接合的
    def update(self, learning_rate):
        pass
    

# 全连接层实现类
class FullConnectedLayer(object):
    def __init__(self, units=-1, activator=None):
        '''
        构造函数
        input_size: 本层输入向量的维度
        output_size: 本层输出向量的维度
        activator: 激活函数
        '''
        #self.input_size = input_dim
        self.output_size = units
        self.activator = activator
        # 权重数组W  
        self.weights=np.ones(1).astype('float32') # / np.sqrt(self.output_size)
        # 偏置项b
        self.weights_bias = np.zeros(1).astype('float32')
        # 输出向量
        #self.output = np.zeros((self.output_size, 1))
        self.initializied=False
    def forward(self, input_array):
        '''
        前向计算
        input_array: 输入向量，维度必须等于input_size
        '''
        # 式2
        #print 'input_array',input_array.shape
        if not self.initializied:
            h,w=input_array.shape
            self.input_size = h
            self.weights = np.random.uniform(-0.1, 0.1,(self.output_size, self.input_size)).astype('float32')
            self.weights_bias = np.zeros((self.output_size, 1)).astype('float32')
            self.initializied=True
        self.input = input_array
        self.output = self.activator.forward(np.dot(self.weights, input_array) + self.weights_bias)           #gnumpy 加速
        
    def backward(self, delta_array):
        '''
        反向计算W和b的梯度
        delta_array: 从上一层传递过来的误差项(也就是本层的误差）(没有乘以激活函数的导数)
        '''
        # 为了得到本层的误差，delta_array 需要乘以本层激活函数的导数
        delta_array=self.activator.backward(self.output)*delta_array
        self.delta = np.dot(self.weights.T, delta_array) # 传到上一层的敏感度向量（每个神经元有一个敏感度） #gnumpy 加速
        self.W_grad = np.dot(delta_array, self.input.T) # 本层权重梯度
        self.b_grad = np.sum(delta_array, axis=1,keepdims=True) # 本层偏置的梯度
      
        
    def update(self, learning_rate):
        '''
        使用梯度下降算法更新权重
        '''
        self.weights -= learning_rate * self.W_grad
        self.weights_bias -= learning_rate * self.b_grad
        
# Sigmoid激活函数类

class SigmoidActivator(object):
    def forward(self, weighted_input):
        return np.tanh(weighted_input)
    def backward(self, output):
        return (1 - np.power(output, 2))
# Relu激活函数类
@jit
def Reluf(input,out):
    if input.ndim==4:
        n,c,h,w=input.shape
        for ni in xrange(n):
            for ci in xrange(c):
                for i in xrange(h):
                    for j in xrange(w):
                        if input[ni,ci,i,j]>0:
                            out[ni,ci,i,j]=input[ni,ci,i,j]
                        else:
                            out[ni,ci,i,j]=0
    if input.ndim==2:
        h,w=input.shape
        for i in xrange(h):
            for j in xrange(w):
                if input[i,j]>0:
                    out[i,j]=input[i,j]
                else:
                    out[i,j]=0        
    #return out

class ReluActivator(object):
    def forward(self, weighted_input):
        return weighted_input*(weighted_input>0).astype(np.float32) #矢量版 效率最高
    def backward(self, output):
        return (output>0).astype(np.float32) #矢量版
# softmax激活函数类
class SoftmaxActivator(object):
    def forward(self, weighted_input):
        tmp=np.exp(weighted_input)
        return tmp/np.sum(tmp, axis=0, keepdims=True) #矢量版
    def backward(self, output):
        return 1 #(output>0).astype(np.float) #矢量版

class IdentityActivator(object):
    def forward(self, weighted_input):
        #return weighted_input
        return  weighted_input
    def backward(self, output):
        return 1
    
    
# 神经网络类
class Sequential(object):
    def __init__(self):
        '''
        构造函数
        '''
        # 所有的层
        self.layers = []
    
    # 用于添加层
    def add(self,layer):
        self.layers.append(layer)
    
    def predict(self, sample):
        '''
        使用神经网络实现预测
        sample: 输入样本
        '''
        output = sample
        for layer in self.layers:
            layer.forward(output)
            output = layer.output
        return output
    def train_epoch(self, labels, data_set, rate,mini_banch):
        '''
        训练函数
        labels: 样本标签
        data_set: 输入样本
        rate: 学习速率
        epoch: 训练轮数
        mini_banch: 小批次
        '''
        batchs=data_set.shape[0]/mini_banch-1
        print 'total samples= %d   batch size= %d  batchs=%d'%(data_set.shape[0],mini_banch,batchs+1)
            # 前 banchs 个banch
        for j in range(batchs):
            t0=time.time()
            self.train_one_sample(labels[:,j*mini_banch:(j+1)*mini_banch], data_set[j*mini_banch:(j+1)*mini_banch], rate)
            #print 'batch ',j,time.time()-t0
            # 最后一个banch
        self.train_one_sample(labels[:,batchs*mini_banch:], data_set[mini_banch*batchs:], rate)
            
    def train_one_sample(self, label, sample, rate):
        self.predict(sample)
        self.calc_gradient(label)
        self.update_weight(rate)
    def calc_gradient(self, label):
        delta = ( self.layers[-1].output-label )
        for layer in self.layers[::-1]:
            layer.backward(delta)
            delta = layer.delta
        return delta
    def update_weight(self, rate):
        for layer in self.layers:
            layer.update(rate)
            
    def save_model(self):
        '''保存模型参数到硬盘文件'''
        f=file("model.npy", "wb")
        for layer in self.layers:
            np.save(f,layer.weights) #保存权值
            np.save(f,layer.weights_bias) #保存偏置
        f.close()
    
    def load_model(self,modelfile='model.npy'):
        '''从硬盘文件加载模型参数'''
        f=file(modelfile, "rb")
        for layer in self.layers:
            layer.weights=np.load(f) #载入权值
            layer.weights_bias=np.load(f) #载入偏置
            layer.initializied=True
        f.close()
    
    # 训练模型
    def train(self,x,y,batch_size=10,epochs=100,rate=0.001,trained=False):
        # 如果已经训练好 直接加载参数
        if trained:
            print 'trained model'
            self.load_model()
            return
        
        # 不然就训练
        print 'training model'
        epoch=0
        loss0=1000
        rate_now=rate
        for i in xrange(epochs):
            t0=time.time()
            epoch += 1
            self.train_epoch(y, x, rate_now,batch_size)
            #print 'epoch %d finished' % ( epoch)
            if epoch % 1 == 0:
                #print network.layers[-1].output.shape
                print 
                print '############################################################'
                print self.layers[-1].output[:,0:4]
                loss=(-np.log(self.layers[-1].output)*y[:,-self.layers[-1].output.shape[1]:]).sum()/self.layers[-1].output.shape[1]
                if loss<loss0:  #误差变小时才保存
                    self.save_model()
                print 'epoch = %d  loss= %f  time= %f '%(epoch,loss,time.time()-t0)
                loss0=loss
                t0=time.time()

    # 计算正确率            
    def evaluate(self, test_data_set, test_labels):
        error = 0
        total = len(test_data_set)
        predict=self.predict(test_data_set).argmax(axis=0)
        for i in range(total):
            if test_labels[predict[i]][i] != 1.0:
                error += 1
        return float(error) / float(total) 